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We study the numerical solution of the nonrelativistic Schrodinger equation for two-electron atoms 
in ground and excited S states using pseudospectral (PS) methods of calculation. The calculation 
achieves convergence rates for the energy, Cauchy error in the wave function, and variance in local 
energy that are exponentially fast for all practical purposes. The method requires three separate 
subdomains to handle the wave function's cusplike behavior near the two-particle coalescences. The 
use of three subdomains is essential to maintaining exponential convergence and is more compu- 
tationally efficient than a single subdomain. A comparison of several different treatments of the 
cusps suggests that the simplest prescription is sufficient. We investigate two alternate methods 
for handling the semi-infinite domain, one which involves a sequence of truncated versions of the 
domain and the other which employs an algebraic mapping of the semi-infinite domain to a finite 
one and imposes no explicit cutoffs on the wave function. The latter prescription proves superior. 
For many purposes it proves unnecessary to handle the three-particle coalescence in a special way. 
The presence of logarithmic terms in the exact solution is expected to limit the convergence to 
being nonexponential but the only clear evidence of that is the rate of convergence of derivatives 
near the three-particle coalescence point. Higher resolution than achieved in this work will ulti- 
mately be needed to see its limiting effect on other measures of error. As developed and applied 
here the PS method has many virtues: no explicit assumptions need be made about the asymptotic 
behavior of the wave function near cusps or at large distances, the local energy {'Htp/tp) is exactly 
equal to the calculated global energy at all collocation points, local errors go down everywhere with 
increasing resolution, the effective basis using Chebyshev polynomials is complete and simple, and 
the method is easily extensible to other bound states. As the number of collocation points grows, 
the method achieves exponential convergence up to the resolution tested. This study serves as a 
proof-of-principle of the method for more general two- and possibly three-electron applications. 

PACS numbers: 31.15.ac,03.65.Gc,02.70.Jn,02.60.Lj,02.30.Jr 



I. INTRODUCTION 

The nonrelativistic, two-clcctron atom (H^, He, Li+) 
is the simplest "hard" problem in quantum mechanics. 
It involves strong electron-electron correlations, nontriv- 
ial symmetry considerations, and single as well as double 
continua. Many different solution techniques have been 
developed and applied over the past 80 years. A thor- 
ough understanding of this simple system is important 
not only because of its direct relevance to experimental 
studies in atomic physics but also because the best meth- 
ods of solution may suggest generalizations applicable to 
multielectron and/or multiatom systems. 

Our own interest in this problem arose from investi- 
gating bound-free and free-free opacity of the negative 
hydrogen ion H^. As first conjectured by Wildt 
gives the greatest contribution to opacity in the atmo- 
sphere of the Sun and many other stars. The photo- 
absorption cross section of H~ is known to an accuracy 
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of a few percentage points [3] but little attention has been 
devoted to H~ in less-than-ideal circumstances (high den- 
sity, high magnetic field, etc.) of relevance to astrophys- 
ical applications. We sought a first-principles approach 
that would allow "exact" calculations of initial and fi- 
nal states as part of these investigations and were led to 
reconsider this classic problem. 

Ideally, there would exist a simple method capable of 
handling any two-electron state in the presence of a nu- 
cleus with any angular momentum whether bound or 
free. In practice, many individual methods have been 
formulated each having somewhat more specific goals. 
A common starting point, for example, is finding the 
ground-state energy for zero total angular momentum. 

It is not possible in a single article, let alone an intro- 
duction, to review the full range of methods that have 
been developed and explored. We can briefiy compare 
the strengths and weaknesses of a few select approaches 
by assessing each in terms of the generality (is it appli- 
cable to all states or just the ground state?), the capa- 
bility of achieving an exact solution of the nonrelativistic 
Schrodinger equation in the limit of infinite nuclear mass 
(is it in principle capable of finding an exact solution [in 
the aforementioned sense] or are there intrinsic approx- 
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imations?), the degree of tuning required (is it straight- 
forward to apply or does it require an enhghtened guess 
for, say, the choice of basis functions?) and, of course, 
the computational effort for a given level of accuracy. 

The asymptotic rate of convergence of some error i?„ 
as a function of the number n of basis functions, grid size, 
etc. is of central importance in evaluating a numerical 
method. To characterize the convergence rate the defini- 
tions of Boyd 3 are used in this article and reproduced 
here. The algebraic index of convergence, k is defined as 
the maximum k so that 

lim \Rn\n^ < oo. (1) 

If k is finite then i?„ converges algebraically. The sim- 
plest example of algebraic convergence is an error i?„ oc 
1/n'^. If k is infinite then i?„ converges exponentially. 
This latter category is subdivided into three cases de- 
fined by the value of 

Z=linil^. (2) 

n— f oo TL 

If / is zero, a finite positive number, or infinite, the rate is 
subgeometric, geometric, or supergeometric, respectively. 
For example, if the error R„ oc exp(— n™), the conditions 
< m < 1, m = 1, and m > 1 correspond to sub- 
geometric, geometric, and supergeometric convergence, 
respectively. 

A. Variational method for two electrons and 
nucleus: ground state 

The first numerical explorations of two-electron ground 
states adopted the approach of minimizing the global en- 
ergy. Once Hylleraas determined that only three coordi- 
nates were needed to represent the wave function for S 
states he carried out such variational calculations (prior 
to the advent of computers) Q. Pekeris and coworkers 
did the first high precision calculations on comput- 
ers, expanding the wave function in terms of Laguerre 
polynomials of linear combinations of the interparticle 
distances times an appropriate exponential falloff. They 
determined the energy of H" to eight decimal places and 
that of He to nine, calculations that were the gold stan- 
dard for several decades. 

Variational methods have been highly successful at 
calculating extremely precise eigenvalues of the ground 
state of two-electron atoms. Indeed, eigenvalue energies 
have been calculated to numerical accuracy — at least 42 
digits — that far exceeds the accuracy of the underlying 
physical description based on nonrelativistic equations 
of motion [9l-[l9||. A clear strength of the general vari- 
ational approach is that intrinsic approximations to the 
Hamiltonian operator need not be made. The principle 
drawbacks are related to the difficulties inherent in the 
selection of the basis: it should be complete so that con- 
vergence to the exact solution is possible and efficient so 



that finite numbers of elements do a good job represent- 
ing the wave function. 

Significant progress in choice of the basis for the two- 
electron problem has taken place. The inclusion of new 
functions (e.g., logarithmic terms) typically motivated by 
known limiting forms of the wave function improves the 
rate of convergence (tI. [qI. [TTI. [20| - |2^ . Furthermore, with- 
out special additions, some bases are simply incapable 
of representing the exact solution [2^ [2^. Klahn and 
Morgan have shown that there are examples where the 
expectation value of an operator (i.e., r'^ with A: > 6) 
converges to the incorrect value or diverges even if the 
basis is complete. In their example the basis cannot ac- 
curately represent the derivatives of the hydrogenic solu- 
tion at r = [131. When employing such bases, one must 
always check that the physical property one is calculating 
is converging properly. 

Schwartz [28| surveyed the convergence rate of the er- 
ror in the ground-state energy eigenvalue achieved by 
many different strategies for basis set selection. His re- 
sults for the error may be expressed as a function of n, 
the total number of basis functions selected according 
to a well-defined procedure. The error generally con- 
verges algebraically with index, 1.5 < k < 8.3 (cx n~'^), 
depending on the basis. The range in k highlights the 
significance that a good choice of basis can have on the 
asymptotic convergence of a calculation. One basis set, 
which included a single power of a logarithm, appeared 
to converge exponentially fast as cr"" with a in the range 
0.51-0.54. Such exponential behavior is often assumed of 
variational methods if the basis can accurately describe 
the behavior of the wave function everywhere. That is, 
the basis includes functions which have the same analytic 
and nonanalytic behavior as the exact solution. 

Loosely speaking, even when convergence is assured, 
the accuracy of the variationally inferred wave function 
(by many different measures) is much less than that of the 
energy eigenvalue. Parts of the wave function that have 
a small effect on the total energy arc not well-constrained 
by lowering the energy. An alternative strategy to min- 
imizing the global energ y is to minimize the variance in 
the local energy instead |29l-l3l| . This approach can pro- 
duce better local values of the wave function but leads 
to nonlinear minimization problems which are more diffi- 
cult to handle numerically (minimization of the variance 
in local energy with respect to parameters in the trial 
wave function) but still tractable because one need not 
calculate the global energy at each step. 



B. Variational method for excited states 

Variational methods d, d, [13, [H, HI] have been suc- 
cessful at calculating precise excitation energies. In gen- 
eral, variational methods extend naturally to excited, 
bound states whenever the variational parameters enter 
in a linear fashion. It is then straightforward to find mul- 
tiple eigenstates of the linear system. The more highly 
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excited the state the less converged the energy is, espe- 
cially if the basis was optimized in order to reproduce the 
features of the ground state only [oij . 

C. Fourier spectral expansion 

Griebel and Hamaekers [131 developed a Fourier ex- 
pansion method for multidimensional quantum mechan- 
ical systems. They apply the hyperbolic cross trunca- 
tion to their basis and show that for smooth solutions 
the exponential convergence rate is not dependent on the 
dimensionality of the system. They calculate the ener- 
gies of several different systems including hydrogen and 
helium. Unfortunately, they fail to achieve exponential 
convergence because the cusps were not properly treated, 
and hence their highest resolution runs for hydrogen and 
helium are only good to about 2 and 10%, respectively. 

D. Specialized methods for two-electron systems 

Haftel and Mandelzweig and later other collaborators 
[35l - l4l| have presented an exact treatment of two-electron 
atoms that begins by factoring out the correct cusp be- 
havior and posing the problem in terms of the remaining 
part of the wave function. This piece which is contin- 
uous up to first derivatives is expanded in terms of hy- 
perspherical harmonics yielding a set of coupled ordinary 
differential equations for coefficients which are functions 
of the hyperradius. The method fully accounts for the 
asymptotic behavior near coalescence points and yields 
results with energies good to one part in 10^. The hy- 
perspherical harmonic expansion converges more quickly 
than if the cusp is not explicitly accommodated for, but 
remains algebraic because of the higher-order discontinu- 
ities. The method accurately determines bound excited 
states as well. 



E. Direct solution of partial differential equation 
for bound and continuum states 

Most of the bound-state techniques mentioned thus 
far are unsuitable for calculating continuum-states. In 
fact, continuum state calculations rely on totally differ- 
ent variational methods. The main ones are i?-matrix 
p^ . Schwinger variational [i^, and the complex Kohn 
variational [4j| methods. Accuracy for these methods 
lags far behind that of bound-state calculations. 

Roughly speaking, the source of some of these diffi- 
culties is related to describing the wave function over 
an infinite volume while simultaneously controlling the 
errors of greatest significance. Typically, linear varia- 
tional methods are equivalent to spectral expansions of 
the wave function. The control one has over the accu- 
racy of an approximate description of the wave function 
is indirect via the choice of the expansion. Instead, one 



may be motivated for both bound and continuum prob- 
lems to consider solving the partial differential equation 
directly on a grid where a greater degree of local control 
is possible. 

Finite difference methods (FDM) jisl - lii l and Finite el- 
ement methods (FEM) |49l454l| represent the solution and 
the differential equation on a discrete grid. The FDM 
grid is usually evenly spaced with derivatives calculated 
to some small (usually second) order. FEM uses subdo- 
mains, concentrating grid points where more accuracy is 
needed. Recent work achieves as many as seven decimal 
places in the energy of the ground state b ut p roduces sur- 
prisingly nonsmooth wave functions (53l [HJ. The rate of 
convergence of these methods is limited by the order of 
the representation of derivatives and is always algebraic 
with some small index dependent on the order used for 
derivatives. 



F. Pseudospectral approach 

Some of the above considerations motivate an inves- 
tigation of the pseudospectral (PS) method. Like FDM 
and FEM methods the PS method represents the wave 
function by values on a discrete grid of points rather 
than by coefficients of a spectral expansion. However, the 
points arc selected in a different manner and the deriva- 
tive order increases with grid resolution. Roots of Ja- 
cobi polynomials are chosen in order to make the asymp- 
totic rate of convergence of an analytic function constant 
across the entire finite nonperiodic domain. Such a choice 
also has the advantage that exponential convergence can 
be lost only by nonanalytic behavior within the domain. 
By contrast, an equispaced grid is sensitive to singular- 
ities nearby in the complex plane and can lead to di- 
vergences when interpolating near the endpoints (Runge 
phenomenon). Of all the Jacobi polynomials, Chebyshev 
polynomials vary the least over [—1,1] and hence produce 
the smallest residual from the PS method. The mathe- 
matical theory of nonsmooth functions is not well devel- 
oped and precise convergence rates are usually calculated 
empirically [ssl ]. 

The PS method [1, [H, [HI] has seen successes in many 
fields including fluid dynamics [s^], relativistic astro- 
physics [11] , and numerical relativity [H, [G^I • When the 
underlying solution is smooth, the PS method typically 
requires less computational run time and less memory 
than FDM and FEM to achieve comparable precision. 
The method has been applied in quantum mechanics to 
solve the full Schrodinger equation for a single electron 
[gH - IgsI ]. In addition, various simplifications of the mul- 
tielectron Schrodinger equation have been treated, in- 
cluding the Hartree-Fock approximation [64l - l69| , M0ller- 
Plesset perturbation theory [t^I , and density functional 
theory [Hi, III- 

To the authors' knowledge, no one has solved the 
full three-dimensional Schrodinger equation for helium- 
like systems (the "exact" problem) using PS methods. 
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This article implements the method, investigates several 
design choices and calculates ground and excited bound 
S states. The convergence rate is used as the metric 
to characterize different grid choices, alternative meth- 
ods for handling regularity conditions and other practical 
considerations needed for an efficient algorithm. No at- 
tempt to reproduce the ultrahigh precision results of vari- 
ational methods is made. The calculation employs a stan- 
dard Chebyshev basis without any specialized tuning. 
The eigenvalue and eigenfunction problems are solved by 
a standard method. All calculations are done on a single 
processor with a speed of 6 GHz and 8 GB of memory. 

The PS method is expected to be supergeometric on 
a finite computational domain if no singularities exist in 
the solution anywhere in the complex plane, geometric 
if singularities are only outside the domain, and alge- 
braic if singularities exist within domain. If the domain 
is infinite or semi-infinite, subgeometric convergence is 
expected when no singularities are in the domain, and 
algebraic convergence is expected otherwise 0] . 



II. SETTING UP THE PROBLEM 

Let Zi and V| be the position vector and the Laplacian 
of the coordinates, respectively, of the ith electron if i is 
1 or 2 and of the nucleus if i is 3. The nonrelativistic 
Schrodinger equation for a heliumlike system is 



where 



n 



2 V m m M 



where 



V 



Z 



Z 



zi - Z3 



Z2 - Z3 



Zl - Z2 



(3) 



(4) 



Z is the nuclear charge, a = 1 unless the electron-electron 
interaction is suppressed (a = 0), m is the mass of the 
electron and M is the mass of the nucleus. The units are 
e = h ~ l/47reo = 1. The Hamiltonian acts on functions 
of nine dimensions, i.e., three coordinate positions for 
each particle. 

The relative and center of mass coordinates are 
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Define the coordinates 

ri 

ri2 



ri = Zl - Z3 

r2 = Z2 - Z3 

m(zi + Z2) + Mzs 



M + 2m 



|ri| 
|r2| 



and rewrite the Hamiltonian 

H — To + Tcm 



r2|, 



V, 
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ri2 



(12) 
(13) 
(14) 
(15) 



/i = mM/{M + m) is the reduced mass of the electron 
and nucleus, and Vx is the gradient operator with respect 
to the vector x. 

In the center-of-mass frame Tcm niay be dropped bring- 
ing to six the number of nontrivial coordinates on which 
the wave function depends. Because m <C A/, the mass 
polarization term Tmp is often ignored or treated pertur- 
batively. While unnecessary for many methods includ- 
ing PS, we use the infinite nuclear mass approximation 
(A/ = cx)) to facilitate comparison with previous results. 
In units with m ~ 1 (atomic units) the Hamiltonian is 



n 



1 



z 

ri 



Z 

7-2 



a 

ri2 



(16) 



Atomic units are used throughout the rest of this article. 

This operator is elliptic. All boundaries in physical 
space require specification of the function or its normal 
derivative or some combination of the two [t^. In the 
ideal problem, the physical boundary is at infinity where 
the wave function must be zero. The existence of the 
Coulomb potential's singular points at ri = 0, r2 = 0, 
and ri2 = introduces complications in any formal and 
practical analysis. Before the exact nature of the Hermi- 
tian Hamiltonian operator and its spectrum was under- 
stood, Kato [tJ showed that discrete eigenstates existed 
for the specific case of helium. In later work Kato [75| 
showed that the wave function must be finite at the sin- 
gular points (which is also true everywhere else) , and that 
the first derivative of the wave function on the domain 
excluding the singular points is bounded. This result al- 
lows discontinuities in the first derivative at the singular 
points, called Kato cusps. Generally, higher derivatives 
are not bounded at the singular points. 

In any numerical treatment of the Hamiltonian opera- 
tor a decision must be made about how to handle the sin- 
gular points. In a formal mathematical sense, quantities 
at the singularities are well defined only in the limit as 
one approaches the singularity. This creates an effective 
inner boundary about such points on which additional 
conditions on the function and its normal derivative may 
be specified. Such conditions are exploited to guarantee 
regularity in the limit that the excised region shrinks to 
a point. This article assumes that it is correct to excise 
such a point, either explicitly or implicitly. 
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III. COORDINATES AND THE HAMILTONIAN spherical coordinates as: 



The hehunihke atom is made of three particles: a nu- 
cleus and two electrons. Six coordinates are required to 
describe relative positions. Three coordinates describe 
the precise shape and size of the triangle with a particle 
at each vertex, and the other three describe the orien- 
tation of that triangle in space (often taken to be Euler 
angles) . The wave function for S states is completely in- 
dependent of the latter three For nonzero angular 
momentum one first expands the wavcfunction in gener- 
alized spherical harmonics of the Euler angles. Only a 
finite number of terms are needed for a given total angu- 
lar momentum and its z component, and the Shrodinger 
equation becomes a finite set of coupled partial differ- 
ential equations for the remaining three variables (e.g., 
Refs. [Zilzl)- 

Two useful sets of coordinates for the triangle are 
{ri, 7'2, 7'i2} and {^"1, ?'2, ^12}, where ri and r2 are the 
proton-electron distances, ri2 is the electron-electron dis- 
tance, and 9i2 is the angle between the vectors pointing 
to the two electrons. Four additional useful sets of coor- 
dinates {p or X, (f>, C} and {p or x, C, B} arc defined by 



n 



n 

C 

V2sinC 
B 



pcos(p 
ps'mcj) 

— cos 6*12 

v/l + Csin20 
cos 

X = . 

1 + p 

The ranges of these variables are given by: 



(17) 
(18) 
(19) 
(20) 

(21) 
(22) 



< ''I , ^2 , p < 00 

|ri - r2| < ri2 < ri -I- 

< 6*12 < TT 

< 0,C < 7^/2 
-1 < x,C,B < 1. 



f2 



(23) 



where 



and 



Tp + p-'^{% + csc^2(j)Tc) + P'^U (25) 
Tp + p-\Tc+csc^2CrB) + P'^U, (26) 



-^app-^9p (27) 



Tc 

Tb 
U 



{l + xf 



{l + xf{A + x) 
4(1 -x) 



2 cot 



-2{{l-C^)dcc~'2Cdc) 



- ( -(9^^ + 2cot2C9(; 



-2 ((1 - B'')dBB - IBOb) 
— Z CSC (j) — Z sec (f> 

Zs/2 Zs/2 



a 



V2sinC (^[BX] (^[-BXV 



(j[x,y\ = + xsn\2y. 



do. (28) 

(29) 

(30) 

(31) 

(32) 
(33) 

(34) 



(35) 



IV. THE SINGULAR POINTS IN THE 
HAMILTONIAN 

PS methods are very sensitive to discontinuous deriva- 
tives of any order. If such discontinuities exist, the 
method loses its exponential convergence and artificial 
oscillations may occur. The wave function has discon- 
tinuities only at the singular points which thus require 
special attention. Myers et al [7^ discuss these singular- 
ities in detail. Here we reproduce some of their discussion 
for completeness. 



The coordinate x maps the semi-infinite domain to a fi- 
nite domain [9^ . This simple choice works well because 
the wave function is exponentially small at large p for 
bound states which are the topic of interest here. 

After integrating over the Euler angles the volume el- 
ements are 

4 / rir2ri2dri dr2 dri2 
4 / r\r2 sin 9i2dri dr2 d9i2 
J p^ sm^ 2(j)dpd(j)dC 
J p^ sm^2CdpdC dB 
2/{l^sin2 2(j)dxd<j>dC 



(24) 



2/ 



sin^ 2Cdx d( dB. 



The Hamiltonian for 5' states can be written in hyper- 



A. Two-particle coalescences 

There exist three lines corresponding to two-particle 
coalescences: two for the proton and each electron at 
= and (p = 7t/2 and one for the two electrons at C = 0. 
Only one of the proton-electron coalescence lines need 
appear in the numerical domain which takes advantage 
of the explicit symmetry of the spatial part of the wave 
function about = 7r/4. 

Kato [79I I analyzed the discontinuity in the derivative 
of a wave function at two particle coalescence points and 
showed that 



dr 



P-tjQtQji'ir = 0), 



(36) 
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where ■tp is the wave function, r is the particle-particle dis- 
tance, ip is the limit of the average value of the wave func- 
tion on a sphere centered at r = as its radius shrinks 
to zero, fXij is the reduced mass of the two particles, and 
Qi and Qj are the charges of the two particles. 

Pack and Byers Brown [s^] extended the analysis to 
show that the wave function could be expanded in terms 
of hydrogenic solutions. 



Im 



•I J 7 



[0, 



1 



(37) 



where I > 0, |m| < Z, aim is an expansion coefficient, 9 
and (f> arc the usual spherical angles giving the orienta- 
tion of the two particles, and Y™' is the usual spherical 
harmonic. 

These results describe the regularity required at the 
Coulomb singularities. There are three practical ap- 
proaches to making sure the solution has the appropriate 
behavior. 

1. Behavioral 

Assume that local solutions to the Schrodinger 
equation that fail to satisfy Eqs. (pG]) and ([37)) are 
not analytic; assume that the expansion (cardinal 
functions) employed in the numerical treatment is 
incapable of representing this nonanalytic behavior. 
Granted these assumptions, all numerical solutions 
will automatically be regular at the point in ques- 
tion According to Boyd 3, in many contexts 
this approach is sufficient. If the solutions that do 
not satisfy Eqs. p6| and ([37]) have only weakly sin- 
gular behavior, the convergence rate may be slow. 

2. Regularity 

Replace the Hamiltonian at the singular points 
with the Kato cusp conditions without otherwise 
altering the domain. The cusp conditions are 



9-0 
c 



4=0 



dip 



Zpi' ( = 2 



f^(C = o) 



(38) 
(39) 
(40) 

(41) 
(42) 
(43) 



where these two sets are mixed and matched while 
preserving the appropriate symmetry or antisym- 
metry. The choice depends on precisely which state 
one wishes to calculate. 



or 



iPicP-- 


= 0) 


= 




1) 


= 


V'(C = 


0) = 


= 0, 



3. Excision 

Excise a small sphere around the singular points 
and impose boundary conditions on its surface that 
yield the correct behavior at the singularity as the 
sphere shrinks. From equation Eq. (j37p 



iP = Y,bip'^'Pi[C] (l- 
i ^ 



l+l 

Zpcp 
l + l 



0[p'<P' 



o 



2 12 



(44) 
(45) 
,(46) 



where <p = tt/2 — cp, Pi is a Legendre polynomial of 
order /, and bi, ci, and di are unknown constants. 
Define 



^Pi[C]dC 



XI = I ^i'Pi[B]dB, 



and write the conditions as 



Z 

TTT 
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a0 + l-^ + ^^7Ti 
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(47) 
(48) 

(49) 
(50) 

(51) 



These conditions become exact as the excised vol- 
ume shrinks to a point. For PS methods the vol- 
ume should be reduced exponentially with increas- 
ing resolution. The changes do not increase the 
computational cost, but may adversely affect the 
condition number of the matrix. 

The difficulty of implementation increases with num- 
ber on the list. 



B. Three particle coalescence 

The potential is also singular when all three particles 
collide (i.e., when the hyperradius, p goes to zero). The 
behavior of the wave function about this point is much 
less well understood than two-particle coalescences and 
cannot be handled in the same way. Instead of simply 
having a discontinuity in the wave function's first deriva- 
tive (the value of which is finite on both sides of the 
singularity), the second derivative grows logarithmically 
near p = 0. Bartlett [2^ was the first to show that a sim- 
ple Frobcnius type expansion in powers of p about p = 
fails at second order on account of the electron-electron 
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interaction. He suggested that logarithmic terms exist in 
the exact solution of helium. Fock [sil, [s^l introduced an 
expansion of the form 



L»/2J 



n— m— 



(52) 



where e„m are two dimensional functions of the hyperan- 
gles {(j), C} or {(, B} determined through the recursive 
relationship 



[n{n + 4) + A]e„m = 2Ven~i,m - 2Een-2,m - 2{n + 2)(m + l)e„.m+i - (m + 1)(to + 2)e„,m+2, (53) 

I 



and A is the two-dimensional Laplacian over the hyper- 
angles. All e„m with n < 2 are known analytically plus 
a few additional terms with higher n |83|-|86|. Morgan 
proved that the series is convergent everywhere [s^l , and 
it has been shown that variational calculations converge 
faster when a single power of a logarithm is included in 
the basis H. 

Again, there are three basic strategies for a numerical 
scheme. 

1. Behavioral 

Do nothing special and rely on the regularity of the 
cardinal functions. This is an imperfect approach 
since the exact wave function has unbounded sec- 
ond derivatives as p — ^ 0. If the basis set can only 
represent regular behavior as discussed in the case 
of two-particle coalescence it will not produce the 
exact solution. However, since the volume element 
scales as p^dp such inexactness may have negligi- 
ble effect on observables calculated from the wave 
function. 

2. Regularity 

Impose a regularity-like condition at the singular 
point. For the ground state (and many other S 
states) , the first-order solution [8l|, [s^] to the Fock 
equations [53] is 



eoo = c 
eio = c 



{~Z{c 



(54) 
(55) 



where c is a constant given by the normalization. 
These solutions imply either 

9pVlp=o = {-Z(cos0 + sin</.) + |(j[C,<?!)]}v(/3 = O), 

(56) 



which is valid for the ground state, or 

V'(p = o) = o. 



(57) 



Note that this regularity-like condition says noth- 
ing about the second derivatives. For the same rea- 
sons as above this method can never give the exact 
solution at p = 0. 



3. Excision 

Excise a small domain with p < pmin, where pmin 
is the cutoff. A boundary condition can be calcu- 
lated on this inner surface by solving equation [53] 
to relate the wave function and normal derivative 
on the surface [o^ . As the resolution increases, 
more terms must be calculated so that the trun- 
cation error in the Fock expansion equals the er- 
ror due to having finite resolution in the numerical 
calculation. The basis set expansion in the bulk re- 
mains completely regular. While exact, the method 
is complicated and will be pursued at a later time. 



C. Infinite separation 

The domain of the ideal problem extends to infinity. 
The bound-state wave functions fall off exponentially. 
The outer boundary condition must be approximated in 
the numerical method. There are several approaches. 

1. Behavioral 

Rely on the regularity of the cardinal functions to 
exclude exponentially growing solutions as p — >■ oo. 
One must map the semi-infinite domain to a finite 
one to use Chebyshev collocation points or work 
with semi-infinite functions like Laguerre polyno- 
mials. 

2. Regularity 

Again map the domain to a finite one but replace 



the Schrodinger equation at p = oo by 
ip{p = oo) = 0. 



(58) 



Note that if one includes the endpoints using 
the behavioral method, this method is effectively 
the same as the behavioral condition because the 
Schrodinger equation reduces to Eip = at a; = — 1 
(p = oo). 

3. Excision 

Excise the region with p > pmax where Pmax is the 
cutoff and impose a suitable boundary condition. 
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As in the three-particle coalescence, one may de- 
velop a more and more accurate representation at 
fixed Pmax and/or an approximate condition at in- 
creasing pmax- It is easiest to set 



or 



Ipip = Pmax) = 



dp 



(59) 



(60) 



and vary Pmax, which is what is done in this article 
when using this method. In a PS numerical scheme 
one should vary pmax oc n.-^^^ for large n where n is 
the radial resolution (see Appendix \Q . 



D. Collinearity {B or C = ±1) 

The coefficients multiplying the second derivatives 
with respect to B and C at B,C = ±1 go to zero. In 
ordinary differential equations this allows irregular so- 
lutions that behave as linear combinations of Legendre 
functions of the second kind. Regularity of the cardinal 
functions excludes such solutions. Since the Schrodinger 
equation at these points contains no infinities it does not 
matter if the grid includes these points. The partial dif- 
ferential equation is parabolic along this boundary. So 
no boundary conditions or regularity conditions need to 
be given. 



V. THE PSEUDOSPECTRAL METHOD 

Boyd 3 , Fornberg [s^ , PfeifFer et al [s^ , and the third 
edition of Numerical Recipes [56j cover the PS method in 
detail. A brief review of some aspects pertinent to our 
work follows. 

The main advantage of this method is that it provides 
exponentially fast convergence for smooth solutions. Un- 
like finite difference and finite element algorithms, all 
derivatives are calculated to higher and higher order with 
increasing resolution. 

It is also noteworthy that the grid points are clustered 
more closely near the boundary of a domain than in its 
center. With this arrangement the representation of a 
function and its derivative is more uniformly accurate 
across the whole domain than is possible using an equal 
number of equidistant points. Finite difference and finite 
clement methods typically use an equal-spaced grid and 
the derivatives are less accurate at the edge than at the 
center. 

Let Ud be the number of coordinate dimensions and 
Ni the resolution in the zth dimension. The differen- 
tial equation is enforced at nt ~ YYHi collocation or 
grid points chosen to be the roots or extrema of a Ja- 
cobi polynomial of order Ni in each dimension. Boyd's 
recommendation that one use Chebyshev polynomials to 



generate the grid points in lieu of special circumstances 
is followed here 

The derivatives in the i"^ direction are calculated to 
Nf^ order in terms of the function values at the collo- 
cation points. To illustrate this it is useful to define the 
cardinal functions: 



N 



1=1 



X — X 

x^ — x^ 



(61) 



where the a;"s are the collocation points and i and j are 
superscripts not exponents. These functions have the 
property that 



Cf[x']^5]. 



(62) 



The PS representation of a function at an arbitrary n^- 
dimensional position {xi, . . . , Xn^) is expanded as 



"0 ~ i^Ni,...,N^^ = 



Ni,...,N„ 



I Xl , . . . , Xyj 



(63) 

(using the Einstein summation convention) in terms of 
its grid values and the cardinal functions 



F, 



^31: 

Ni,...Mr, 



i^Nu...,N,Jx{\...,x'nf] (64) 



(65) 



Such an expansion is equivalent (up to an exponentially 
small error for smooth functions) to a spectral one. 



(66) 



where 



Gj = X{n,Ax^] (67) 
T/iGj J|w[a;,]da;, (68) 

rid 

« ^^^'•••■'".Gj[xl\...,x:r/]n^'»., (69) 

fe=i 

j means ji , . . . , , Uj- arc orthonormal polynomials cho- 
sen here to be Chebyshev polynomials, w is the weight 
function over which they are orthonormal, and Vi^ are the 
corresponding quadrature weights. Note, the collocation 
points are also quadrature points which allow exponen- 
tial convergence of the quadrature. Derivatives of ip are 
approximated by differentiating Eq. (j63p . To this end, it 
is useful to introduce differentiation matrices. 



{DN)]^d.,.Cf[x%=, 



(70) 



This method is readily applied to linear eigenvalue 
problems as arise from the Schrodinger equation [96j 



[n - E)i> = 0, 



(71) 
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where H is the Hamiltonian operator, tp is the wave func- 
tion, and E is the energy eigenvalue. Discretizing on the 
grid gives the tensor equation 



H 



where 



H 



31,32, 



31,32, 



d^3 



^F3l,. 



32, 



\ ^1 *2 



■ J •'-lid 



Unlike finite difference methods, the tensor H 



dense. Write this as H, 



k{ii,i2,...,i„.) 



31,32,-- 

or for short Hj' , 



(72) 



(73) 



where 



l{3l,32,---,jn 

k and I are one-to-one functions mapping the set of rid 
indices to the lowest nt positive integers. This recasts the 
tensor as a large matrix so that standard matrix methods 
can be employed. 

One way to carry out the mapping employs the Kro- 
necker product as follows. If H is given by 



n 



/n,...,^„,[2;l,...,a;„J(9:,J'l •••(a:,„J*"^ 



(74) 

where is a function coefficient, then the matrix. 

El is given by 



— ^ ^ fii,...,i 




(75) 



where x'^ is the Ti^-dimensional vector of coordinates with 
indices that map to k, Nj is the number of grid points 
in the j*^ direction, ij is an exponent, and Dn^ is the 
differential matrix based on Nj points. 



VI. PSEUDOSPECTRAL CONVERGENCE FOR 
NONSMOOTH FUNCTIONS 

To make appropriate design algorithmic choices it is 
important to investigate how the PS method handles 
nonsmooth behavior of solutions. This section explores 
the convergence of truncated cardinal function expan- 
sions to cusps and logarithmic terms and then employs 
a toy model that illustrates how the triple coalescence is 
expected to influence the numerical results. 



A. Kato cusps 

Consider the ground state of the hydrogen atom with 
wave function 



(76) 



In Cartesian coordinates, there is a discontinuity in the 
first derivative at the origin. 



dip 

2:^0+ OX 



y=z=0 



, y dip 
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X 

V- 

10 12 14 



FIG. 1: (Color online). The logarithm base 10 of ^RMsbi] 
(blue circles) and ^RMslSa] (red crosses) with solid blue and 
dashed red fits, respectively. See Appendix iBl for fitting func- 
tions. 



In spherical coordinates no discontinuity exists for r > 
0. All the derivatives at r = are well defined and a 
PS code has no problem exponentially converging toward 
the correct answer. The essence of this observation can 
be seen by considering the one-dimensional exponential 
functions 



9im 

52 N 



^(x+l) 



(78) 
(79) 



on the domain — 1 < a; < 1 with weight (analogous to 
the three dimensional hydrogen atom). As a measure of 
error between the function / and its cardinal expansion 
truncated at order n define 



RMS 



^/^:^^(^/N-ECfN/M) dx. (80) 



Figure[T]compares (5rms [ffi] t° '^rms [52] as a function of n. 
Evidently the cusp is poorly represented compared to the 
smooth function at a given n. The PS representation of 
the cusp converges algebraically while the representation 
of the smooth function converges supergeometrically (see 
Appendix IB] for fits). 

The basic strategy in more complicated problems is to 
adopt a coordinate system with a radial-like coordinate 
at each cusp. For two-electron atoms no global coordi- 
nate system exists with the desired property at each of 
the three separate two-particle coalescences. This article 
uses three individual but overlapping domains to guar- 
antee appropriate treatment near each coalescence point. 



y=z=0 



B. Logarithmic terms 

Consider the one-dimensional function 
(77) /[x] = (^1 + \og[p[x]]^ e-^W, (81) 
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FIG. 2: (Color online). The logarithm base 10 of the error in 
f[x] using cardinal functions at p = 10~® (green pluses), p = 
lO"'^ (red crosses), and p = 10* (blue circles) with solid blue, 
dashed red, and dotted green fits, respectively. See Appendix 
|B]for fitting functions and method. 



where 



1 - X 



p[x\ = 



1 



(82) 



Here p £ [0,oo), a; G [—1,1] and / are analogous to the 
hyperspherical radius, its algebraic transformation and 
the heliumlike wavefunction respectively. As in the 
full three-dimensional problem, the presence or absence 
of the logarithmic terms is controlled by a, which can be 
set to or 1. 

There are two types of errors considered here: inter- 
polation error and operator error. These are different 
sorts of error, but qualitative features (e.g., exponential 
or algebraic convergence) are expected to be the same. 



1. Interpolation error 

The pointwise error between / and its truncated ex- 
pansion is 



N 



(83) 



i=l 



where refers to the ith grid point. 

Figure [2] shows that the behavior of A/ at three dif- 
ferent values of p. For each value the apparent rate of 
convergence starts out exponential before becoming alge- 
braic at large N. The algebraic convergence known with 
the highest accuracy in Fig. [2] is for p = 10~^, which 
asymptotically goes as i/7V^-^^='=°-°^ (see Appendix IB|) . 
This algebraic behavior is expected when trying to rep- 
resent a nonanalytic function (log[/9]) with an analytic 
basis. Such behavior disappears if a is set to zero. 

The onset of algebraic convergence varies from N « 
AO to N 80 as p, moving away from the singularity, 
increases by 10 orders of magnitude. The error at the 



transition is < 10~^. Typically, energy errors vary as 
the square of wave function errors and would already 
be very small compared to relativistic corrections. This 
calculation shows that it is possible to get precise values 
with apparent exponential convergence before reaching 
the asymptotic algebraic regime. As a practical matter, 
one may never reach the latter limit. 



2. Operator error 

In order to estimate the error in the eigenvalue, it 
would help to have a one-dimensional toy eigenvalue 
problem with an eigenfunction similar to the function in 
Eq. ([5T|) . It is impossible to construct a one-dimensional 
eigenvalue problem with solutions that have logarithmic 
singularities without explicitly introducing such singular- 
ities into the differential operator. So here a more limited 
test problem is used. Instead of solving for an eigenvalue, 
the error in the operator is measured. This would con- 
tribute to the eigenvalue error along with the error in the 
wave function. 

Let AH be the difference between the true Hamiltonian 
and the Hamiltonian constructed from PS differentiation 
matrices. The associated energy error is (tplAT-Llip) . The 
aim is to construct an analog of the integrand of the 
energy error and use it to assess pointwise and integral 
errors. 

Construct a differential operator V similar to the full 
Hamiltonian H [see Eq. ((25)) ] but in terms of the coordi- 
nate X, 



where 



P2[x] = 

Pi[x] = 
Po[x] = 



_{i + xY 

8 

(1 + .t)3 f A + x 
1 

p[x] ■ 



l-x 



(84) 

(85) 
(86) 
(87) 



Note, the first two terms are identical to the operator Tp 
and the last term is a Coulomb potential. 
The corresponding matrix operator is 

(rfAr)} =P2[a^1(i?Ar)U^w)^? +Pl[a^1pAr)} +Po[x1<5;, 

(88) 

where Einstein's summation convention is used and Sj is 
the Kronecker delta function. The pointwise error on the 
grid and its maximum are 



(AH 



n) = w[x 



{Vf)[x^]-Y.idN))f[x=]] (89) 



AW^-^ = max|(Ai7jv)i, 



(90) 
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FIG. 3: (Color online). The logarithm base 10 of the maxi- 
mum error of a pseudospectral matrix. The dark blue circles 
are for a — 1 and the light red crosses for a = with solid 
blue and dashed red fits, respectively. See Appendix [B] for 
fitting functions and method. 



metric because the calculation is done on a semi-infinite 
domain). The log term does spoil the method's exponen- 
tial convergence. Assuming that the effect is greatest at 
i = i* , the maximum error is dominated by the log term 
when N is greater than about 200 and the error is very 
small. This is exponential convergence "for all practical 
purposes." 

C. Conclusion 

For the interpolation and operator errors, the logarith- 
mic term does not slow convergence unless one is at high 
resolution or interested in small values of p. For those 
cases, one would need to apply the excision method about 
the triple coalescence point in order to retain exponential 
convergence. For most applications the level of precision 
needed is obtained before the algebraic behavior becomes 
apparent. 




20 40 60 80 100 



FIG. 4: (Color online). The logarithm base 10 of the error of 
a pseudospectral matrix near the singularity. The dark blue 
circles are for q = 1 and the light red pluses for q = with 
solid blue and dashed red fits, respectively. See Appendix iBl 
for fitting functions and method. 

where 

w[x] = ^j—^M (91) 
[1+ x)' 

The factor (1 — x)^ /{\ -\- xY comes from the Jacobian. 

Figure |3] shows AHf}'^^, the maximum error anywhere 
on the grid, as a function of n for a = and a ~ 1. 
The decrease appears to be exponential, not unantici- 
pated when a — but perhaps a surprise for a = 1. The 
slopes of the two curves are roughly the same and the 
offset is due to the variation of the magnitude of / with 
a. An explanation is immediately suggested by Fig. |3] 
which shows the error at the grid point closest to the sin- 
gularity {AHnY (here i* refers to that point). The data 
for a = 1 is well fit by an algebraic rate of convergence 
(■j^^^io.36±o.08) at large N while a = has an approx- 
imately exponential fall-off (the convergence is subgeo- 



VII. NUMERICAL DOMAINS AND 
COLLOCATION POINTS 

Because the two electrons are identical particles the 
full wave function is antisymmetric. In the ground state, 
the spins are antisymmetric and the spatial part is sym- 
metric. The spatial domain may be taken as (/) < 7r/4 
and B > 0. The complete numerical domain is 

-Do: -l<x <1 0<(/)<f -1<C<1 

or (92) 
Do:-l<a;<10<C<f 0<S<1. 

A single domain does not allow the proper treatment of 
the two-particle coalescences. Therefore, introduce three 
subdomains to cover Dq using two different sets of vari- 
ables: 

Di: -l<a;<l O<0<i -1<C<1 

D2: ~l<x<l \<(t><l -|<C<1 (93) 

Dz: -l<x<l Q<Q<\ 0<B<1. 

For calculations done on a finite domain, the condition 
— 1 < a; < 1 is replaced by < p < pmax. Cross sections 
of these domains at fixed p are shown in Fig. [51 An 
electron-proton singularity lies in Di, while the electron- 
electron singularity lies in D^,. The radial- like coordi- 
nates in Di {(f)) and {() accommodate the cusps just 
like the usual radial coordinate does in the hydrogen 
atom. D2 fills in the remaining volume. All three do- 
mains have boundaries that touch the triple coalescence 
point (not pictured). 

Consideration of the electron-electron singularity 
shows why the single domain Dq is inadequate. Byers 
Brown and White [88[ showed that the wavefunction can 
be expanded in powers of ri2 about ri2 = 0. Using such 
a coordinate accurately treats the cusp away from the 
triple coalescence point. The expansion in powers of ^ 
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Di and D2 and {A^^, -^s, -^ci = i^*^! '^i in domain D3. 
The total number of grid points is rij = 3 x (2n xnxn) = 
6n^ points. Twice as many points were used in the x 
dimension, an arbitrary choice but one motivated by the 
semi-infinite range of the hyperspherical coordinate and 
by the wave function's logarithmic dependence on the 
hyperspherical radius near the triple coalescence point. 



VIII. BOUNDARY CONDITIONS 



A. Internal boundary conditions 



FIG. 5: (Color online). This is the arrangement of grid points 
of the three domains at a constant value of p in and C co- 
ordinates. Note that the point density becomes larger at the 
boundary of each subdomain and that no grid points sit on 
the Coulomb singularities. The blue circles, red crosses, and 
green pluses belong to domains Di , D2 , and D3 , respectively. 
Di and D2 are rectangular domains, while D3 has the curved 
boundary in cp, C coordinates but is rectangular in B co- 
ordinates. The electron-proton singularity occurs on the left 
side (solid line ai <j) — 0). The entire line corresponds to one 
physical point. The electron-electron singularity occurs at the 
lower right-hand corner (solid disk at0 = 7r/4, C = — 1). A 
line of symmetry falls on the right side (dashed line at tf) — n/A 
where ri = r2). 



is very similar [see Eq 
y2sinC 



(j46p ] or equivalently powers of 
C sm2(j) [q^]. Re-expanding in (f> and 



C coordinates gives derivatives of \/l + C sin 20 with re- 
spect to C and 0, terms that are either infinite or un- 
defined at = 7r/4 and C = -1. This is why the PS 
method fails to converge rapidly using Dq alone. 

Within each domain, grid points are set as follows. 
Let the ith dimension extend from Xi^min to Xi^max and 
have Ni collocation points. These points are the roots 
or antinodes plus endpoints of the N!^^ order Chebyshev 
polynomial |98j stretched to fit length Axi = Xi^max — 
2^1, min- The jth poiut (for j = 1, 2, . . . , Ni) of dimension 
i is 



Ax,; 



where 



Vi 



[N, - J + A)7r 



(94) 



(95) 



and A is or 1/2 for nodes or antinodes plus endpoints, 
respectively |99[. In this article, nodes are generally used 
except when explicit boundary conditions are needed at 
both endpoints, Xi^min and x^^max. 

Potentially each dimension and domain could have its 
own Ni but in this paper the x direction is set to be twice 
as large as the other two dimensions and all are varied in 
lockstep. That is, {N^, Nc, N^} = {2n, n, n} in domains 



It is necessary to ensure continuity of the wave func- 
tion and its normal derivative at internal boundaries. 
There are two ways in which the subdomains can touch: 
they can overlap or they can barely touch. For clarity, 
consider a one-dimensional problem with two domains. 
Let the first domain be domain 1 and the second be do- 
main 2 with extrema Xi,,nin < a;2,min < Xi^^ax < X2,max, 

where the 1 and 2 now refer to domain number. The 
first case corresponds to X2,min < xi^max and the second 



to X2. 



xi. 



For both cases, exactly two 



conditions are need to make the wave function and its 
derivative continuous. The simplest choice for the first 



V'l[a^l,inax] 
■0l[a^2,min] 

and for the second case is 



'02[xi,max] 
-02 [X2,min], 



(96) 
(97) 



-02 [a; 
d d 



-02 [x* 



(99) 



dx dx 
For multidimensional grids, the situation is analogous. 
The conditions are applied on surfaces and the deriva- 
tives are normal derivatives at the surface. On a discrete 
grid, a finite number of conditions are given which, in the 
limit of an infinitely fine mesh, would cover the entire 
surface. There is a great deal of freedom in the selec- 
tion of the points but in this article the edge of a domain 
has one constant coordinate so there is a natural choice. 
Conditions are imposed at the points of the finite mesh 
formed by varying all the other coordinates (in general, 
these are not collocation points). In other words, the 
matching points lie at the intersection of the coordinate 
lines normal to the surface with the surface itself. The 
positions of the crosses in Figs. [S] and [7] illustrate where 
the matching occurs when the domains overlap and when 
they just touch. 

For touching domains, the black and white crosses in 
Fig. [7] are used. Note that four (three) crosses are de- 
fined by the coordinate lines in D2 {Di). At the set of 
four crosses, function values are equated, and at the set of 
three crosses, normal derivatives are equated. In general, 
function values (derivatives) are equated at points stem- 
ming from the subdomain with a greater (lesser) density 
of points along the boundary. 
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FIG. 6: The intersection (gray) of domains D\ (white) with 
black grid points and Da (black) with white grid points. The 
boundary points are depicted as black and white crosses and 
are connected via black and white lines to the grid points that 
they replace. 



D, 




FIG. 7: Barely touching domains D\ (white) with black grid 
points and D2 (black) with white grid points. The bound- 
ary points are depicted as black and white crosses and are 
connected via gray lines to the grid points that they replace. 

For overlapping domains, the function values are 
equated at the black and white crosses in Fig. [S] which 
lie on two separate surfaces. The points are selected in a 
manner similar to that of touching domains, i.e. in terms 
of the intersection of the coordinate lines in Di and D3 
with the surface. 

In all cases, values and derivatives at all points are 
calculated using Eq. (j63p and are ultimately linear com- 
binations of the grid point values. 

B. Symmetry and regularity conditions 

For this problem there are two types of boundary con- 
ditions on the boundary of the numerical domain: the 
symmetry condition from electron exchange and the reg- 
ularity conditions imposed near singular points. 

The symmetry condition is related to the total spin of 
the two electrons, S. If 5 = (5* = 1) the wave function 
is symmetric (antisymmetric) about (/) = 7r/4ori3 = 



and the normal derivative (value) of the wave function is 
equal to zero. This condition is enforced at all the points 
on the boundary that have the same p and C coordinates 
as grid points in D2 or the same p and C coordinates in 
D3. This gives 2 x 2n x n = An? boundary conditions. 

Regularity conditions are imposed as boundary condi- 
tions at four two-dimensional surfaces: p = 0, p = 00, 
(/) = 0, and C = 0. These are similar in form to the 
symmetry condition except involve linear combinations 
of derivative and value. Depending on which type of 
conditions are given at singular points (behavioral versus 
regularity or excision) are used, there are to lOn^ condi- 
tions. These conditions replace an equal number of equa- 
tions. The particular equation replaced is the one that 
stems from enforcement of the discrctized Schrodingcr 
equation at the collocation point nearest to the bound- 
ary at which the condition applies. 

The most complicated type of boundary condition 
arises when a region is excised about a two-particle coa- 
lescence. First, one must project out terms proportional 
to each Legendre polynomial by performing an integral 
over C or B. This can be done by quadrature over the 
grid points in those dimensions. For example, Eq. (|47p 
turns into 

ii[p\<j>i] = (100) 

k 

where Wk are the quadrature weights. Then Eq. (|49|) 
becomes for j = 1 (the excision boundary) 

= {<l>'{DN,)i + [-1 + ^) <5^) ii[p\<i^% (101) 
which is NpNc conditions (0 < / < Nq — 1 and 1 < i < 

C. Incorporating boundary conditions into the 
matrix problem 

All of the above boundary conditions are expressed as 
a linear combination of the function values at the grid 
points equal zero. In matrix form 

where ipi {1P2) is a vector of the rib (rii) wave function 
values at all the boundary (interior) points, the boundary 
condition matrix has been broken into an ni, by Uf, matrix 
Bi and an rib by i^i matrix B2, and Ub + n.i = rit. For 
the case where an endpoint is not a collocation point, the 
grid point nearest to the boundary, at which an explicit 
boundary condition is given, is considered as a boundary 
point. All the points near where behavioral boundary 
conditions are given are not included in this definition. 
These points are the ones that give rise to the first Ub 
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rows in Eqs. (|102[) and (|103|) . Note, that this ordering 
was chosen for clarity in this section. 

There is also the Hamiltonian matrix equation 

^ « ■' 

where the Hamiltonian matrix has also been divided into 
four matrices: Hn, H12, H21 and H22- 

So there are rit + rib equations and unknowns (ipi 
and "02) as well as the eigenvalue. One could approxi- 
mately solve these equations with singular value decom- 
position [5^, but it is much faster to simply discard the 
first rib rows of the Hamiltonian matrix and incorporate 
the boundary conditions into the remaining eigenvalue 
problem: 

{H22 - H2iB^^B2 - E)^2 = 0, (104) 

where Bi has an inverse because all of its rows are lin- 
early independent (otherwise more than one boundary 
condition would have been specified for a given bound- 
ary point). Calculating the inverse is not too computa- 
tionally expensive because rib ^ nt- Then solve for ipi 
afterwards with 

V-i = -B^^B2^2- (105) 

IX. MATRIX METHODS 

In one dimension the Hamiltonian matrix for the PS 
method is dense but in three dimensions with three do- 
mains the number of nonzero elements scales as 24n* 
out of a possible 36n^. The boundary condition matrix 
is also sparse with 8ri^ non-zero elements out of 48ri^ 
(for the simplest case where behavioral conditions are 
used whenever possible). Therefore, any attempt to solve 
these equations should take advantage of these memory 
savings. 

Equation (|104p is solved by the method of inverse it- 
eration [56j after shifting the eigenvalue with an approx- 
imately known value. In cases where the exact eigen- 
value is not known a priori, one solves the full eigenvalue 
problem for a low resolution case first and then at each 
successive iteration shifts the eigenvalue using the result 
of the previous iteration. Because the Hamiltonian ma- 
trix is not symmetric, a complex eigenvalue may occur. 
There is no theoretical reason prohibiting the numerical 
eigenvalue from containing an imaginary part at finite 
resolution but, in fact, none were generated for n > 5. 
Of course, the imaginary part contributes to the error 
which must converge to zero. 

The above solution method can yield highly oscillatory 
wave functions which appear to diverge on the boundaries 
of the computational domain. These nonphysical wave 
functions do not satisfy the first rib rows of Eq. (|103p 
and arise as an artifact of solving a subset of equations 



TABLE I: A list of the different cases that are compared in 
this section. Exc refers to the first excited Estate. Grd refers 
to the ground state. Nd is the number of domains. B, R, and 
E refer to behavioral, regularity, and excision, respectively. 
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Grd 


3 
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of the overdetermined system. They are easily identified 
and rejected and in no way affect the true solution. 

The entire calculation for n = 14 took only about 20 
minutes on a 6-GHz machine. Memory needed to solve 
the linear equations was the limiting factor because in- 
verting the equation has requirements scaling as . The 
generalized minimal residual (GMRES) algorithm [s^ 
might reduce the memory requirements of the solution 
of the linear equations that arise in the inverse iteration. 

For simplicity of coding, the above calculations were 
done using Mathematica [90|. Care was taken to use 
predefined functions whenever such choice was more effi- 
cient. 



X. RESULTS 

This article is an exploration of the PS method as ap- 
plied to heliumlike systems, not an attempt to improve 
the energy eigenvalues for bound states. That has al- 
ready been done to a higher precision than will ever be 
needed The focus here is on showing the PS 

method works in a new application and assessing its con- 
vergence properties. 

Tabic H] gives a list of runs used in this section to 
discuss the effects of the Coulomb terms, energy level 
(ground or excited), computational domains and numer- 
ical methodology on the convergence of the solution to 
the two-electron problem. 



A. Convergence in energy 

In this section, the energy error means the difference 
between the numerical energy eigenvalue at finite resolu- 
tion and the exact energy eigenvalue of the nonrelativis- 
tic infinitc-mass-nuclcus Hamiltonian. When no analytic 
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FIG. 8: (Color online). The convergence of the energy of H~, 
case A (red pluses); two non-interacting electrons in the field 
of a proton, case B (green stars); the ground state of He, case 
C (blue crosses); and its first excited S-state, case D (black 
circles) as a function of grid resolution n, with dotted red, dot- 
dashed green, dashed blue, and solid black fits, respectively. 
See Appendix [Bl for fitting functions and method. 



FIG. 9: (Color online). The convergence of the energy of, 
case E, H~ using only one computational domain with the 
electron-electron interaction on (blue circles); case F, one do- 
main with the interaction off (red crosses); and case A, three 
domains with the interaction on (green pluses) with dashed 
blue, dotted red, and solid green fits, respectively. See Ap- 
pendix |B] for fitting functions and method. 



value exists, highly precise variationally calculated values 
are used [9l-[l9|. 

The energy errors for H~, H~ with the electron- 
electron interaction turned off, the ground states of He- 
lium, and the first excited 5 State of helium (cases A,B,C, 
and D) are shown in Fig. [S] The first important result 
is that the energies appear to converge in an approxi- 
mate exponential fashion. Since these are not variational 
calculations there is no reason to expect monotonically 
decreasing energy errors. Detailed inspection of the solu- 
tions suggests that the kinks in the graphs are discrete- 
ness effects. That is, the precise positioning of the grid 
points has a large effect on the magnitude of the error. 

A potentially significant issue is the impact on con- 
vergence of logarithmic terms present in the Fock expan- 
sion. Prior authors have been able to calculate very accu- 
rate energies without adverse effects of the infinite second 
derivative at p = 0, but the PS method is sensitive to sin- 
gularities anywhere within the domain. Figure |H] includes 
a calculation of the H~ system with the electron-electron 
interaction turned off, altering the exact solution and re- 
moving the logarithmic term. The rate of convergence 
is comparable in all cases suggesting that the influence 
of the logarithmic term on convergence is subdominant 
for n < 14. This conclusion agrees with the analysis in 
subsection IVIBI 

Ideally, the numerical method should handle states 
other than the ground state. Figure |8] shows the con- 
vergence of energies for the ground state and the first 
excited S state of helium. The important result is that 
the convergence of both calculations is approximately ex- 
ponential with a similar rate. 

The relative sizes of the magnitude of the error at fixed 
grid size for H~, He and excited He are roughly consis- 
tent with the general expectation set by the difficulty in 



resolving the solution's small-scale structure. Errors for 
ground state He are larger than H~ because the expo- 
nential length scale for falloff of the He wave function is 
smaller than that of H~ ; errors for the excited state of 
He are larger than the ground state of He because the os- 
cillatory length scale of the excited state is smaller than 
the exponential length scale of the ground state. 

Figure [5] shows the impact on convergence of us- 
ing a single numerical subdomain, Dq, versus three, 
{Di, D2, D3}. The single domain had one third as many 
points as the computation with three domains. However, 
the resolution in the x direction dominates the conver- 
gence and in that dimension the resolution is identical. 
Domain Dq has radial-like coordinates near the electron- 
proton cusp but not near the electron-electron cusp. One 
anticipates slower convergence in the energy using Dq. 
Comparison shows that two interacting electrons on three 
domains (case A) or two noninteracting electrons on Dq 
(case F) have similar exponential rates of convergence. 
On the other hand, two interacting electrons on Dq (case 
E) converge more slowly. This result shows the multiple 
grids are essential for achieving superior convergence, and 
that the electron-electron cusp drives this requirement. If 
the single domain data were adjusted to account for the 
fact that they used fewer points, they would be shifted 
to the left by a factor of 3^^^ « 1.44. This would not 
affect the conclusion that using three domains is more 
efficient because the single domain solution is only alge- 
braically convergent starting at about n = 8. Using three 
subdomains is more efficient because the work involved 
in the calculation has the same scaling whether one or 
three subdomains is used. 

Figure [TO] presents a comparison of calculations hav- 
ing the full semi-infinite domain (case A) to those with 
a finite cutoff in p (case H). The scaling of the cutoff 
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FIG. 10: (Color online). The convergence of the energy of, 
case H, H~ doing the calculation on a finite domain (blue 
circles); and, case A, semi-infinite domain (red pluses) with 
dashed blue and solid red fits, respectively. See Appendix iBl 
for fitting functions and method. 

Pmax ^/n imposed in case H is derived in Appendix |A] 
by balancing the error due to finite resolution from the 
numerical scheme with errors introduced by truncating 
the bound state. The figure shows that the semi-infinite 
calculation fairs better. This is a consequence of the two 
different sets of assumptions used to distribute the points. 
The grid points in the semi-infinite scheme are more often 
found where the wave function is large. Half the points 
have p < 1 (a; > 0) because is the center of the x dimen- 
sion. By comparison, half the points have p < Pmax/2 m 
the finite calculation. The number of points where the 
wave function is large is smaller in this latter scheme. Al- 
though the semi-infinite strategy is more effective, noth- 
ing can be said about the optimal strategy because other 
distribution methods were not considered. The main ad- 
vantage of the method is simplicity since there are no 
adjustable parameters. 

Figure [Tl] presents a comparison of the different ways 
of handling the regularity of the wave function at the 
two-particle coalescence points. The simplest method, 
relying on the regularity of the Chebyshev polynomials 
(case A), does as well or better than the other methods 
(cases I and J). 

Figure [T^] compares two ways of handling the wave 
function at p = 0, case A, relying on the regularity of the 
Chebyshev polynomials (behavioral) and case G, directly 
specifying a logarithmic derivative (regularity). The lat- 
ter method is slightly better but both have roughly the 
same convergence rate. 

B. Convergence in local energy 

Another useful measure of convergence is the local en- 
ergy, 

Eioc - (106) 



FIG. 11: (Color online). The convergence of the energy of 
H~ using three different methods of ensuring regularity at 
the two-particle coalescence points: case A, relying on the 
regularity of the Chebyshev polynomials (green pluses); case 
I, using the Kato cusp condition as a regularity condition (red 
crosses); and case J, excising the singularity (blue circles) with 
green dotted, solid red, and dashed blue fits, respectively. See 
Appendix iBl for fitting functions. 
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FIG. 12: (Color online). The convergence of the energy of 
H~ using two different methods of ensuring regularity at the 
three-particle coalescence point: case A, relying on the regu- 
larity of the Chebyshev polynomials (blue circles); and case 
G, using the Fock condition to specify a logarithmic derivative 
(red pluses) with dashed blue and solid red fits, respectively. 
See Appendix iBl for fitting functions and method. 

which is constant only for an exact eigenfunction ?/; of 
Hamiltonian "H. Throughout this subsection all analysis 
and data refers to case A. 

The difference between the local energy and the nu- 
merically evaluated eigenvalue E gives a local measure of 
the error in i\j in a particular calculation. Define 

ASioc = Sloe - E. (107) 

For the PS method, Aii^ioc is zero at all grid points (sub- 
ject to limits of finite precision arithmetic). Nonzero 
differences exist between grid points. Figure [T^] illus- 
trates the convergence of local energy at four different 
points. Of the four points, the error in local energy is 
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FIG. 13: (Color online). The convergence of the local en- 
ergy of H~ at a four points in the domain: the center of the 
computational domain (black circles), near the triple coales- 
cence point (blue crosses), near the proton-electron coales- 
cence point (red pluses), and at large p (green stars). Their 
geometric fits are given by the solid black, dashed blue, dotted 
red, and dot-dashed green lines, respectively. See Appendix 
|B]for fitting functions and method. 
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FIG. 14: (Color online). The error in the local energy of H~ as 
a function of x with ri — 2r2, and C = — 1 (the electrons are 
on the same side of the nucleus) at four different resolutions: 
n — 5 green dot-dashed, n = 8 red dotted, n = 11 blue 
dashed, n = 14 black solid. 



lowest at the point in the center of the computational 
domain {{p,(f>,C} = {l,n/8,0}). It is larger in mag- 
nitude near the singularities ({p, 0, C} ~ {10~^, 7r/8, 0} 
and {1,10^^,0}) because near these points Eioc -> oo 
for any nonexact -0. However, the geometric fits show 
that the rate of convergence is approximately the same 
at all three of those points. A different behavior is seen 
at {p, 0, C} — {15, 7r/8, 0}. The error is roughly constant 
over much of the graph, but begins to decrease at high 
resolution. This is not surprising given the fact that there 
are only a few grid points at such large hyperradius. 

Figure [T3] displays the convergence of the local energy 
as a function of x at fixed angular coordinates. For a 
perfect exponential decrease in local energy error, the 
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FIG. 15: (Color online). The error in local energy of 
plotted at different values of the angle Q\2 with r\ = 2r2 and 
at resolution n — and at p = 1 (solid black), p = 0.1 
(dashed blue), p = 0.01 (dotted red), p — 0.001 (solid greed), 
and p — 0.0001 (dot-dashed purple). 



curves would be equidistant from each other. This is 
approximately true throughout the domain except near 
X = ±1 (small and large p). 

If the numerical solution is considered to be trust- 
worthy where the local energy error is less than some 
threshold (e.g. |A£:ioc| < lO'^), then the wave func- 
tion is well represented in an intermediate range of p 
(10~^ < p < 10^ '^) but not near the triple coalescence 
point nor at infinity [lOOj ]. 

At large p {x w —1), the true wave function falls off 
exponentially (in fact, with respect to the a; coordinate 
it falls off even faster). The PS method represents the 
exponential in terms of a polynomial. When one extrapo- 
lates using the polynomial to a; = — 1, the wave function 
is small but nonzero (the exact value should be zero). 
However, the Hamiltonian acting on the polynomial is 
guaranteed to be zero because every coefficient in the 
Hamiltonian operator has a factor of (1 -t- a;). So 



lim ~ 0, 



(108) 



and Ai^ioc ^E, a constant at large p as seen in Fig. 
UM Detailed inspection of the data near x = —1 suggests 
that for any finite p there exists a resolution above which 
the solution becomes trustworthy. 

The local energy behavior near the triple-coalescence 
point {p ~ 0) is of special interest as a probe of the 
wave function's nonanalytic behavior. Figure [15] displays 
Ai?ioc as a function of ^12, the angle between the two 
electrons, for fixed r2/ri and for a number of choices of 
p, following a similar figure from Myers et al In 
this small p regime, the terms that dominate the Hamil- 
tonian are the kinetic energies in the various directions. 
Each of these, individually, scales as l/p^ [see Eq. p5|) ]. 
For the exact solution, these terms cancel each other, 
but for almost any solution which is not exact, the lo- 
cal energy scales as 1/p^. This scaling is shown in Figs. 
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FIG. 16: (Color online). The root-mean-square average 
Cauchy error is plotted with increasing resolution for three 
cases: H~ with noninteracting electrons (blue circles); H~ 
with interacting electrons (red crosses); and helium with in- 
teracting electrons (green pluses) with dashed blue, solid red, 
and dotted green fits, respectively. See Appendix [Bl for fitting 
functions and method. 



FIG. 17: (Color online). Pointwise differences in the wave 
function evaluated at p = for increasing resolution for three 
cases: H~ with noninteracting electrons (blue circles); H~ 
with interacting electrons (red crosses); and helium with in- 
teracting electrons (green pluses) with dashed blue, solid red, 
and dotted green fits, respectively. See Appendix[B]for fitting 
functions and method. 



[U and [151 Again detailed inspection of the data near 
X ^ \ suggests that for any finite nonzero p there exists 
a resolution above which the solution becomes trustwor- 
thy. Furthermore, Fig. [T3l shows that there is no sign of 
the convergence rate being slowed due to the logarithmic 
terms at this resolution. 



C. Cauchy errors 

Throughout this subsection all data refers to cases A, 
B, and C. The Caueliy error is a measure of the differ- 
ence between numerical solutions with different resolu- 
tion. One such measure is the normed quantity 



The true satisfies 
1 = 



(109) 



(110) 



but integrating ijjn all the way to p = oo would diverge. 
This is a consequence of having small but nonzero errors 
at /9 = oo in the value of il^n- An upper limit p = 10 
is adopted in the normalization of VAi and calculation of 
A„. It is arbitrary but encompasses most of the physical 
extent of the solution. The Cauchy error in any subin- 
terval of the full interval must converge. To the extent 
that the error in the interval calculated is dominant, the 
rate of convergence can be assessed. 

Figure [16] gives A„ as a function of resolution while 
Fig. [ITj gives the pointwise difference at p = where 
the wave function is maximum. Both plots show that 
convergence is approximately exponential. 




FIG. 18: (Color online). The root-mean-square error in the 
logarithmic derivative evaluated at p = with increasing reso- 
lution for three cases: H~ with noninteracting electrons (blue 
circles); H~ with interacting electrons (red crosses); and he- 
lium with interacting electrons (green pluses) with dashed 
blue, solid red, and dotted green fits, respectively. See Ap- 
pendix |B| for fitting functions and method. 



D. The logarithmic derivative at the triple 
coalescence point 

Throughout this subsection all data refers to cases A, 
B, and C. The only direct evidence that the convergence 
of the solutions is slowed by the logarithmic terms in 
the exact solution comes from evaluating the logarithmic 
derivative with respect to p at p = 0. The exact value is 



= -\ (-^(cos( 



(HI) 



The root-mean-squarc error is 
where 

j df7 = ^ d(l)sm^2(j) J dC. (113) 

Figure [18] displays (5r,ms- Turning off the electron- 
electron interaction, the convergence is noticeably faster. 
The important conclusion is that the wave functions' 
convergence is indistinguishable from exponentially fast, 
while the convergence of its derivatives, at least near 
/9 = is slower. Of course, the second derivative with 
respect to p is infinitely wrong at p = 0. It converges to 
a finite value, and the exact value is infinite. 

XI. CONCLUSION 

This article demonstrates the application of PS meth- 
ods for solving the nonrelativistic Schrodinger equation 
for a system with two electrons. The method successfully 
handled both ground and excited S states of heliumlikc 
systems. 

The rate of convergence for most properties measured 
was indistinguishable from being exponentially fast. Lo- 
cal errors decrease in the same manner. 

The choice of variables in the vicinity of the two- 
particle coalescence and the use of multiple, overlapping 
domains are the critical requirements. These are impor- 
tant so the PS method can represent the analytic form 
of the solution near all the two-particle cusps and en- 
sure a more efficient algorithm. In other respects the 
most straightforward choices work well. For example, 
grid points are determined by the roots of Chebyshev 
polynomials, which experience shows generally produce 
the best convergence in PS methods [s, [Hj . Behavioral 
boundary conditions (no explicit regularity conditions) 
are sufficient to handle the wave function in the vicinity 
of the coalescence points and also produce convergence 
as good as or better than the other possibilities tested. 

The energy eigenvalue found by the variational method 
converges most efficiently when basis functions which be- 
have like the exact solution are included but this selec- 
tion process can be time-consuming and problematic. Of 
course, much higher precision than reported here was ob- 
tained long ago by variational methods. The PS method 
has the advantage in new and possibly also in more com- 
plex applications of not needing the same sort of special- 
ized tuning that has benefited variational calculations. 
Although this article does not attempt to reproduce the 
ultra-high-precision results achieved by variational meth- 
ods it strongly suggests that the PS method will ulti- 
mately prove to be a superior approach for reaching such 
results in systems with a small number of electrons. 
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FIG. 19: Energy error from truncation as a function of p, 
with a fit to the points from pmax = 15 to pmax ~ 25. 



The local energy is not directly controlled when total 
energy is minimized. Local energy minimization schemes 
exist and have the advantage that excited states are 
found at local minima of the variance in local energy in- 
stead of just the ground state as is true for the standard 
variational method [3l| . However, they lead to nonlinear 
problems which may be difficult to handle numerically 
(minimization of the variance in local energy with re- 
spect to parameters in the trial wave function), but still 
tractable because one need not calculate the energy at 
each step. By contrast, the PS method controls local en- 
ergy while the numerical solution remains a linear one. 
At the same time, the PS method is superior in terms 
of its convergence rate to other direct partial differential 
equation solvers (grid based methods) such as finite dif- 
ferencing and finite element which also control the local 
error. 

We plan to extend the method to calculate non-S states 
and continuum two-electron states to compute the pho- 
toabsorption bound-free cross sections with both initial 
and final states evaluated with the same methodology. 



Appendix A: Setting pmax for the Truncated Domain 

When truncating the domain at a finite p, it is wise 
to balance the error produced by the cutoff with that 
given by the finite resolution. The former was studied by 
calculating the energy of case H as a function of pmax- 
Figure [19] shows that difference between the truncated 
energy and the correct value of the energy falls off as 
an exponential. The calculation was done with n = 14. 
We see at large values of Pmax that this finite resolution 
ruins the exponential behavior. There is a minimum at 
Pmax = 31. This is probably where the resolution er- 
ror happens to cancel the truncation error. At larger 
Pmax, the resolution error dominates. Therefore, only 
the points with 15 < Pmax < 25 were used for the fit, 
logj^o \AE\ = Apmax + B. A and B were found to be 
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Fit: log, JAE| = -0.6387 n + 1.1344 




FIG. 20: Energy error from finite resolution as a function of 
resolution at pmax ~ 20 with a fit. 



-0.2089 and -0.4313, respectively. 

In order to measure the effect of finite resolution pmax 
was fixed at 20 and the difference between the energy at 
resolution n and 14 was plotted in Fig. [501 The er- 
ror from resolution effects should increase when pmax 
is increased because the density of points goes down. 
So the resolution error is assumed to be of the form 



logio \^E\ = Cn/p 

n 



D. C/20 and D were found 



to be —0.6387 and 1.1344, respectively. 

Setting the two errors equal to each other yields the 
formula 



Pn 



= -3.7476 + 7.8100Vn + 0.2297. 



(Al) 



It is technically possible to get better energies as was 
shown in Fig. [19] (due to error cancellations) , but taking 
advantage of that kind of effect is fine-tuning. 



converging to zero as n goes to infinity, consider fitting 
functions of the forms 

/gcomM = 01(02)" (Bl) 
UM - ai^"' (B2) 

/sup/subM = 01(02)""^ (B3) 
/geom/algM = Oi (a2 )" "I- a3n°^ (B4) 

which are geometric, algebraic, supergeometric (03 > 1) 
or subgeometric (03 < 1), and mixed geometric and al- 
gebraic fits, respectively. A fitting method is used: 



Qn 



f[r 



(B5) 



where Q„ is the quantity Q evaluated at resolution n and 
the sum over n goes from the even numbers from 8 to 100 
for the one-dimensional models, 5 to 14 for the Cauchy 
errors, and from 4 to 14 for all others, is minimized 
with respect to all for the fit which is most reasonable 
on theoretical grounds. Errors in the a,; are estimated by 
calculating 



X 



(B6) 



at the minimum of . Of course large error in the values 
of amplitudes ai and 03 for the mixed geometric and al- 
gebraic fit imply that the other parameters that multiply 
that amplitude may be meaningless. 
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Appendix B: Fits 

Numerical rates for convergence are summarized in Ta- 
ble [III The method is as follows: let Q„ be the quantity 
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TABLE II: The convergence fits of quantities, Q, displayed in figures throughout the article. 
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